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Abstract 

We propose a new evolution equation for the gluon density relevant 
for the region of small xb- It generalizes the GLR equation and allows 
deeper penetration in dense parton systems than the GLR equation 
does. Lhis generalization consists of taking shadowing effects more 
comprehensively into account by including multigluon correlations, 
and allowing for an arbitrary initial gluon distribution in a hadron. 
We solve the new equation for fixed a s . We find that the effects of 
multigluon correlations on the deep-inelastic structure function are 
small. 



laenenOsuryall . cern. ch 
2 levin@laf ex . cbpf .br ; levin@ccsg.tau .ac.il 

3 On leave from Theory Department, St. Petersburg Nuclear Physics Institute, 188350, 
St. Petersburg, Gatchina, Russia 



1 



1 Introduction 



Our objective in this paper is to derive a new evolution equation describ- 
ing the behavior of large partonic densities. The fact that parton densities 
increase as the Bjorken scaling variable xb decreases follows directly from 
linear evolutions equations such as the GLAP |J or BFKL equation, and 
has been experimentally observed by both HERA collaborations 0. They 
observe a powerlike behavior of the deep inelastic structure function 



with ujq ~ 0.3 — 0.5. This behavior is predicted by the BFKL equation and is 
not inconsistent with the GLAP equation. (Such powerlike growth at small 
x can be imitated by a solution of the GLAP equation with a distribution 
that is constant in x at a low initial scale Qo-) 

From a physical point of view such behavior is inconsistent with unitarity. 
This fact necessitates a change in the evolution equation in the region of small 
xb- The first attempt, more than ten years old, to write down a new equation 
led to the nonlinear GLR equation Q. See || for a more extensive review. 

We will explain in this paper that the GLR equation only includes two- 
gluon correlations in the parton cascade. Multigluon correlations should be 
essential to solving the small xb problem, at least theoretically @, 0- R is 
the aim of this paper to take such correlations into account. 

In section 2 we derive and explain the limitations of the GLR equation. 
We consider multigluon correlations in section 3 by employing the relation 
between these correlations and high twist contributions to the deep-inelastic 
structure function. In section 4 we suggest the evolution equation which takes 
these correlations into account. Some particular solutions and a comparison 
with the GLR equation are discussed in section 5, while we discuss the general 
solution of the new evolution equation in section 6. We summarize and 
conclude in section 7. 



In the region of small Xb and large Q 2 we face a system of partons at mutually 
small distances (thus the QCD coupling a s is still small), but dense enough 
that the usual perturbative QCD (pQCD) methods cannot be applied. The 




(1) 



2 The GLR Equation 
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physics that governs this region is nonperturbative, but of a different nature 
than the one associated with large distances. The latter corresponds to the 
confinement region, and is usually analyzed using lattice field theory or QCD 
sum rules. In contrast, we encounter here a situation where new methods 
must be devised to analyze such a dense relativistic system of partons in a 
non-equilibrium state. We need, in fact, new quantum statistical methods 
to describe the behavior of such a system and to chart this unknown region. 
We are unfortunately only at the beginning of this exploration. 

On the upside, we can approach this kinematical region in theory from 
the pQCD region, and assume that in a transition region between pQCD and 
high density QCD (hdQCD) we can study such a dense system in some detail. 
To illustrate what new physics one might anticipate in this transition region 
let us look at deep-inelastic scattering and compare with the pQCD results 
for this process. We can expect the following phenomena in the transition 
region: 

(i) As xb decreases the total cross section o"(7*iV) grows and, near the 
border with the hdQCD domain, becomes even comparable with the geomet- 
rical size of the nucleon o"(7*iV) — > a em irR%. Here the cross section should 
be a smooth function of InQ 2 . 

(ii) Although the parton language can be used to discuss the main prop- 
erties of the process, interactions between partons become important, es- 
pecially the annihilation process. This interaction induces screening (a.k.a. 
shadowing) corrections. 

(iii) In this particular kinematical region such screening corrections are 
fortunately under theoretical control. They modify however the pQCD linear 
evolution equation. The correct evolution equation now becomes nonlinear. 

First some nomenclature. In this paper we will only deal with the gluon 
density in the nucleon, xbG(xb,Q 2 ), which can be measured fairly directly 
at small xb in e.g experiments on diffraction dissociation or heavy quark pro- 
duction in deeply inelastic processes or even directly from the deep-inelastic 
structure function F 2 (xb,Q 2 ) ||- We will thus also refer to it as the gluon 
structure function. 

We will now give a simple derivation of this nonlinear equation, based 
on physical concepts. As we stated, the main new processes that we must 
consider at high density are parton-parton interactions. To incorporate these 
in our physical descriptions we must identify a new small parameter that 
controls the accuracy of calculations involving these interactions. Such a 
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parameter is 

a 



Q 

where p is the gluon density in the transverse plane 



W = 7^-P> ( 2 ) 



x B G(x B ,Q 2 ) , oN 

o = ^ ' (3) 

Here characterizes the area of a hadron which the gluons populate. It 
is the correlation length of gluons inside a hadron. Naturally, this radius 
must be smaller than the radius of a hadron (proton). Since this paper is 
mostly devoted to the discussion of purely theoretical questions we will not 
specify further the value of R^. However it should be stressed that Rn is 
nonperturbative in nature: all physics that occurs at distance scales larger 
than i?7v is nonperturbative. 

The first factor in (0) is the cross section for gluon absorption by a par- 
ton in the hadron. Hence W has the simple physical meaning of being the 
probability of parton (mainly gluon) recombination in the parton cascade. 
The unitarity constraint mentioned in the introduction can be represented 
as § 

W < 1 . (4) 

Thus W is indeed the small parameter sought. The parton cascade can be 
expressed as a perturbation expansion in this parameter. This perturbative 
series can in fact be resummed |4j] and the result understood by considering 
the structure of the QCD cascade in a fast hadron. 

There are two elementary processes in the cascade that impact on the 
number of partons. 

splitting (1 — > 2); probability oc a s p ; (5) 

annihilation (2 — > 1); probability oc a 2 s d 2 p 2 oc a 2 p 2 , 

where d is the size of the parton produced in the annihilation process. In the 
case of deep- inelastic scatterinq d 2 ~ l/Q 2 - 

When xb is not too small only the splitting of one parton into two counts 
because p is small. However as xb — ► annihilation comes into play as p 
increases. 
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This simple picture allows us to write an equation for the change in the 
parton density in a 'phase space' cell of volume A In A In Q 2 : 

d 2 p a s N c af77T 2 

p t^tP > ( 6 ) 



dln-^dlnQ 2 vr Q 2 

X B ^< 



where iV c is the number of colors. In terms of the gluon structure function 

d 2 x B G(x Bl Q 2 ) a s N c 2 2 2 

-x B G{x B ,Q ) - , v ; „ 2 {x B G[x B , Q )) . (7) 



din ^ din Q 2 7T a v "'^ 7 Q 2 i? : 



This is the so-called GLR equation [[|. To determine the value of 7 and 
the understand the kinematical range of validity of (0) this simple physical 
description does not suffice; rather one must analyze the process carefully in 
P^-perturbation theory [|] . The result for 7 was found to be 

81 

7 = — for iV c = 3. 
lb 

We would like to emphasize that the main assumption in the above derivation 
was that 

P {2) - P\ (8) 

where denotes the probability for two gluons in the parton cascade to 
have the same fraction of energy x and tranverse momentum (characterized 
by r ~ InQ 2 ). 

By assuming @ we neglect all correlations between the two gluons other 
than the fact that they are distributed in the hadron disc of radius Rn- 
In the large x B region this assumption is plausible because the correlations 
are power suppressed, and the densities involved are small. In the small x B 
region we cannot justify even when it holds for large x B . 

3 Induced multigluon correlations 

In [^J it was shown that the problem is oversimplified when one assumes that 
the probability of annihilation is simply proportional to p 2 in deriving the 
GLR equation. It was found that 



P 2 



p(2) 1 /ttgfiSain^iaJ^ 



5 



where Qq is the initial virtuality in the parton cascade. This ratio increases 
with decreasing x B - Consequently we must take dynamical correlations into 
account, which could change the GLR equation crucially. 

The key to the calculation of parton correlations was suggested by Ellis, 



Furmanski and Petronzio in []10| , and was developed further in [|IT|]. It was 
shown that gluonic correlations are directly related to higher twist contri- 
butions in the Wilson Operator Product Expansion (OPE) arising from so 
called quasi partonic operators. According to the OPE, the gluon structure 
function can be written as 

x B G(x B ,Q 2 ) = x B G^(x B ,Q 2 ) + q^x 2 b G^(x b .Q 2 ) 

...+ ... {QRN ) 2{n . 1} x B G^(x B ,Q 2 )... (10) 

where the n'th term results from the twist 2n quasi-partonic operator. The 
probability density to find n gluons in the cascade with the same x and 
Q 2 can be directly expressed through the n'th term in the above expansion 
(P (B) =x B G^(x B ,Q 2 )/(nR 2 N ) n , P® = p). 

We recently determined the anomalous dimensions 72« of these high twist 
operators [[F] to next-to-leading order in the number of colors iV c . This was 
done by reducing the complicated problem of gluon-gluon interactions to 
rescattering of colorless gluon 'ladders' (Pomerons) in the t-channel. In [[J it 
was shown that this approach works for the case of the anomalous dimension 
of the twist 4 operator. The fact that there is no Pomeron 'creation' or 
'absorption' in the t-channel means that we are dealing with a quantum 
mechanical problem (not a field theoretical one) in the calculation of the 
72n anomalous dimension. The problem then amounts to calculating the 
ground state energy of an n-particle system with an attractive interaction 
given by a four-particle contact term (see Fig.l) of strength A. Its value can 
be calculated to be 

a N 

X = 4^5, (11) 

7T 

where 5 = 1/{N 2 — 1) if one only takes color singlet ladders into account. 
Including the other color states renormalizes 5 to 0.098 ||. This effective 
theory is two-dimensional (the two dimensions corresponding to ln(l/x) and 
ln(Q 2 )) and is known as the Nonlinear Schrodinger Equation. It is well known 
that this model is exactly solvable. Translating the ground state energy of 
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this model into the value of the anomalous dimension led to ^ 



72n = ^-{1 + \(n 2 -l)}, (12) 

where a s = a s N c /ii, to — N — 1, N being the Mellin-conjugate variable to 

x B 

f(N) = / dx B x N B ~ l f{x B ) or f{u) = / dye-^[x B f(x B )}, (13) 
Jo Jo 

where / is an arbitrary function and y = ln(l/x B ). 

The answer (O) is only reliable when 5 2 n 2 /3 <C 1 M- Thus to check 
selfconsistency we must first generalize the GLR equation based on ([TJ) 
and understand what value of n is important for the deep-inelastic structure 
function. If the answer is inconsistent with the condition 5 2 n 2 /3 < 1 we 
must try and find the expression for the anomalous dimension valid for any 
n. 

Let us note that if we neglect the term proportional S 2 in (|12"D we can 
estimate x B G^ n \x B ^Q 2 ) via the inverse Mellin transform 

x n B G {n \x B ,Q 2 ) = — f due (ujy+ ^ )r) M^ n \u,Q 2 = Ql) , (14) 

where y = \n(l/x B ), r = \n(Q 2 /Ql), Q being the initial virtuality in the 
parton cascade. The contour is to the right of all singularities in M as well 
as to the right of the saddle point (u)s) which is given by 

-J-{ U V + l2nr}\u=u s = 0. (15) 
Thus in the saddle point approximation for 5 = 

x n B G^{x B ,Q 2 ) ~ [x B G(x B ,Q 2 )} n , (16) 

in particular = p 2 . We also assume here the factorization of the ma- 
trix element M (n) = (M (1) ) n ; this expresses the physical assumption that 
there are no correlations between gluons other than the fact that they are 
distributed in a disc of radius R^. We will comment more about this as- 
sumption further on. In this sense the GLR equation is only the lowest order 
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approximation in S 2 even if we assume the factorization of the matrix element 
and thus, strictly speaking, is valid only in the limit of a large number of 
colors. The term proportional to 5 2 in ( |T!?D is clearly responsible for induced 
gluon correlations and we take it seriously in this paper. In next section we 
will therefore generalize the GLR equation. 



4 A New Evolution Equation 

The first step in such a generalization is to make an ansatz for P^ using 
the same approach as for the GLR equation, viz. the competition of two 
processes in the parton cascade. Thus, in analogy to the derivation of (BI) we 
write 

d 2 P^(x B ,Q 2 ) (n) 2 a 2 sl ir {n+x)( 2 , . 

dln^dlnQ* ~ C 2n -P (x B ,Q) ~ n-—P (x B ,Q ) , (17) 

X B ^ ^ 

where C*2 n = 72n^- The factor n in front of the second term on the right 
hand side of flTT| ) reflects the fact that in the Born approximation n + 1 gluons 
annihilate in n gluons through the subprocess gluon + gluon — > gluon, which 
corresponds to the two — ladder — > one — ladder transition with the strength 
of the triple ladder vertex 7. There are n such possibilities due to the time 
ordering of gluon emission. Note that the infinitely recursive set of equations 
can be cut off at any level by imposing e.g. P^ = p( n - m )pM. The GLR 
equation is the case P^ = (P*- 1 -*) 2 . Since we operate under the assumption 
that high twists are essential for small enough i^we must however consider 
the whole series in (|T0D, which we now do using eqs. (|17j). 
Let us introduce the generating function 

00 

gixBiQ 2 ^) = £e"V n) , (18) 

n=l 

where = x B G <yn \xB,Q 2 )- Comparing with fllCf ) we see that for the full 
structure function 

x B G(x b , Q 2 ) = Q 2 R 2 N g(x B , Q 2 ,v=~ HQ 2 R 2 N ))- (19) 
The recursive set of equations (|T^) can be summarized in one equation for g: 
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To solve this linear, 4th order partial differential equation in three variables, 
we must impose some boundary and initial conditions, on the Q 2 and xb 
behavior respectively. 

The boundary condition is straightforward 

For T), ln(— ) fixed; g(x B ,Q 2 ,v) lnQ -?°° e v g L LA(x B , Q 2 ), (21) 
x B 

where Qlla is the solution of the standard GLAP evolution equation for the 
leading twist gluon density in leading In Q 2 approximation. 

The initial condition is much more difficult, because we need g(xs = 
x° B ,Q 2 ,r)) for solving (20), whereas experimentally we can only measure the 
structure function, which is at fixed r]. In other words, we need detailed 
information on the gluon distribution in a hadron at large Xb- We can make 
the following suggestion (although others are possible) 

g(x%Q 2 , V )= (22) 

E^-^r-[9LLA(x B ,Q 2 )r = l-exp(-e^4(4,Q 2 ))- 
This can be recognized as the usual eikonal approximation for the virtual 



gluon-hadron interaction | 12| . The virtues of this expression lie in the fact 
that it is simple, and that it has the transparent physical meaning of reflecting 
the assumption that there are no correlations between gluons with x ~ 1 
other than that they are distributed in a hadron disc of radius R^. 

If one replaces the hadron with a nucleus, one can prove such an approach, 
which corresponds to the so-called Glauber Theory of shadowing corrections. 
In the deep-inelastic scattering case an expression of this type was discussed 



by A. Mueller in |2| 



To summarize this section, we have proposed a new evolution equation 
which has two new features over and above the GLR equation: 

(i) It includes induced multigluon correlations. 

(ii) It allows an arbitrary initial condition not necessarily an eikonal one, 
unlike the case of the GLR equation. We recall that the GLR equation has 
been proven only under the assumption of the factorization property of the 
matrix elements, which corresponds to an eikonal initial condition H. By 
including all twists in our evolution equation we overcome the need for a 
reductive initial condition, such as (^). One could e.g. try to solve (p0|) 
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using an initial condition with correlated gluons at large x and study the 
consequences of its evolution, with or without 5. Because such initial corre- 
lations must be very small, and because our main interest lies in comparing 
with the GLR equation, we do not pursue this line of inquiry here. 



5 Approximate Solutions 

In the next section we will discuss the general solution to eq. (|17D . Here we 
will give approximate solutions for various special cases. In particular we 
try to answer the following questions: how does the nonlinear GLR equation 
follow from our present linear equation; what value of n is typically relevant 
in the sum (|18j) and how do the corrections to the ordinary GLAP evolution 
due to (|i~7l) differ from those due to the GLR equation? 



5.1 The GLR Equation from the Generalized Equation 

The first question that arises is how the nonlinear GLR equation is con- 
tained in the linear equation ( |TTD if we neglect the term proportional to 5 2 . 
Specifically, we would like to establish the equivalence of 

to the nonlinear GLR equation. Let us parametrize the solution in the form 

g(x B , Q 2 ,V) = H e v F(y = - In x B , r = HQ 2 R 2 N ) )) . (24) 

Because we dropped the gluon-correlation term from (|l7|) we can impose the 
'no-correlation' initial condition at Xb ~ 1- This is the only initial condition 
the GLR equation allows 0. In the set of "fan" diagrams (Fig.3) that the 
GLR equation sums the initial distribution has the form 



gix^Q 2 ^) = J2e^(-ir-[g LLA (x%Q 2 )r , (25) 



71=1 



which gives 
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The absence of the 1/n! in the above equation compared with (|22| ) is ex- 
plained as follows. The 1/n! in (|22"D enforces the correct time ordering of the 
produced partons in the parton cascade, related to diagrams of production of 
n parton shadows (see Fig. 4 which shows a case in which three parton shad- 
ows are produced). In the fan diagram of Fig.3 we do not have to enforce 
the correct time ordering because it is already included via the vertex 7 Q. 
Thus the initial condition of eq. (|26D just corresponds to the sum of "fan" 
diagrams, of which Fig.3 is the lowest order example, with the assumption 
that there are no correlations between gluons inside the proton. We will 
see in section 6 that the main properties of the full solution of eq. ([H]) do 
not depend on the form of the initial condition. The reduction to the GLR 
equation does however. 

With (|26|) it is now straightforward to check that eq. (|23|) reduces to 
the GLR equation to first order in e v (recall that finally we must put 77 = 
- \n(Q 2 R 2 N ) < 1), with F(y, r) = xG(x, Q 2 ). 



5.2 Relevant Twists in the Solution 

Here we try to follow the recipe mentioned below eq. (|13|): to see if our 
approach is consistent we must determine what values of n are relevant in 
the solution to (^). Should those values not be consistent with 5 2 n 2 /3 <C 1 
then we must generalize the expression (|1^) for 7 2n such that it is valid for 
all n. If they are consistent we can maintain (|i~2|) and thus (f20|). 

However, to determine the relevant n's we need the general solution to 
eq. (^Oj), which is presented in section 6. Here we will perform some rough 
estimates. 

Let us first take a 'worst case' scenario by letting the term proportional 
to 5 2 dominate, and by dropping the first term and last term on the RHS of 
©: 

d 2 g(x B ,Q 2 ,r]) aj^_ d^g_ _ dg_ 

dln^dlnQ 2 3 % 4 drf' 1 ' 

To solve this equation we perform a Laplace transform in the variable rj 

~g(x B ,Q 2 ,p)= f° d^e^gixs.Q 2 ^). (28) 

4 The strength of the three Pomeron vertex 7 was calculated in ref. || using the AGK 
cutting rules 113], which are equivalent to time-ordering. 
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Then 

d 2 g(y,r,p) a s 5 2 2 2 

— tt^ = ~5--p(p -l)wr,p), 29 

ayar 3 

where y = ln(l/xfi), r = 1ii(Q 2 /Qq). The solution to ( p9|) is the Bessel 
function Io(2^J~^p-p 2 (p 2 — l)yr). The most general solution is then 



9(y, r, 77) = / c ^ c-^0(p) J ( 2^ V - l)yr ) , (30) 

where the contour C runs to the right of all singularities in p. The function 
<p(p) is fixed by imposing an initial condition (e.g. ([22])) 

9(0,r, V ) = ^-J c dpe~ p y(p). (31) 

One may write the solution in fact directly in terms of the initial condition: 

g(y,r,v) = / dr}G(y,r,ri-r})g(0,r,ri), (32) 

with the Green's function 



G(y,r,rj-ri') = ^ f dp e^-^I ( 2i/-^p 2 (p 2 - l)yr ) . (33) 



From eq. ( pT8| ) we see that large typical n corresponds to small typical 77, 
which in turn by (|3CiD corresponds to large typical p ("po")- Thus we must 
find 5 2 p 2 < 1 to trust eq. ([TJ). 

Let us investigate (|32|) to find the most relevant values of 77'. The func- 
tion g(0,r,T]') falls down monotonously for 77' — > — 00 (the first term in (|18|) 
dominates) where it behaves as exp(r]'). Next, we can show that the func- 
tion G(y, r,r] — has a maximum at 77 = 77' whose width is of the order 

of (8J 9 ^f-yr) ( - 1 / 2 \ To see this, evaluate (Bl) using the asymptotic form for 



Io(z) ~ e z /^/2ttz(1 + ...). Eq. (|33| ) then becomes (neglecting the unimpor- 
tant non-exponential prefactor) in the asymptotic form 



n( l\ 1 f, -P( V - V ) + 2 1/ ^fV (2 _! )~ 

G(j/, r,rj-rj) = — J dp e V 3 



(34) 
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The equation that determines the saddle point ps here is 
or 

1 7] — 7] 



, n 2(2p|-l) /c^ 2 

- fa - V) + / \Hr- • v r = ( 35 ) 



VpI- 1 + -T^— = - P ' • (36) 



In the dangerous case that ps is large this leads to 



With (34) this becomes 



PS = I- ,2 ■ ( 3T ) 



_l /2 - iiltf 



G(y,r,r]-r}') = [S^^-yrj -e 8 V^-. (38) 

Thus, in summary, either or both the regions rj ~ and rf ~ 77 give dominant 
contributions in (p3|). For 77' ~ we have 



//(//. /•. 77) - < " Xi 
while for 77' < — ' 77 

g(y,r,rj) oc e~ H . 

The first of these can dominate the second only for a restricted range in 77, 
namely 

la 5 2 

\v\<^^--yr, (39) 

which implies, with ( |3"TD 

< 2 . (40) 
Outside of the range (|39f) the typical value of (77 — 7/) is at large 77 of order 
2^1 . yr) (1/2) , yielding 

p s 1 < 1 for 7/r > 1. (41) 

(2y^ • yr)(V2) 
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We conclude, based on the simplified model equation ( p3| ) that typical values 
of p, and thus n, are small, and therefore we should be able to use eq. (|T2|) . 
We will see that this conclusion holds when we consider the full solution in 
section 6. 



5.3 Estimates for Corrections to the GLR Equation 

In this subsection we want to estimate the possible size of corrections to 
the GLR equation. Recall that the GLR equation sums the contributions 
of fan diagrams (Fig. 3) while the generalized equation includes more general 
graphs, which are all of the type shown in Fig.5 at large Q 2 . To estimate the 
correction let us parametrize the full solution as 

g(x B , Q 2 , rj) = $( e"F(y, r) ) + Ag(y, r, 77) , (42) 

where the first term is the solution to the GLR equation found in subsection 
5.1, and the second term is considered to be a small perturbation. Equation 
becomes then for Ag 



d 2 Ag(y,r,r]) d 2 Ag a s 5 2 . d 4 $ d 2 $ . ^„ r rt rt ^ 2l2 „ 
' /; = a ± -i - — ( ) — <f> " \ p ' p' — a F ]e v 

dydr s d V 2 + 3 1 dif ' [ y r s 1 ' 

(43) 

where $' denotes d<&/dt. We neglected the contribution of the dAg/drj — Ag 
term in eq. fl4~3|) since the coefficient in front of this term contains an extra 
power of as which we treat as a small parameter. The initial conditions 
for Ag(y,r,r)) are Ag(0,r, r])— Ag(y,0,rj) = 0, because we assume that all 
boundary conditions have been fulfilled by $. The last term in eq. (|43D can 
be neglected for two reasons: (i) at r\ — — \ia(Q 2 R 2 v ) this term is suppressed 
and (ii) the difference in brackets is small in both the semicassical and the 
EKL M approach. Substitution of the explicit form of $ in (EBT) yields 



d 2 Ag(y,r,r)) _ (FAg aj?_ 12t 2 (t - 1) 

dydr ' ' as drj 2 3 ' (1 + t) 5 ' 1 ' 

where t = e^Fiy^r). We can simplify the equation if we keep in mind that 
the value of rj is large and negative in the deep-inelastic structure function. 
Thus t <1 and 

9^ g (yr, v) = a< ^_ 12 g^ FVv ( 
oyor or] 1 3 
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Now the T) depedence of Ag is trivial: Ag = e 2v AF(y,r). Thus 

^^=4^,0-12^. (46, 

Using similar techniques as in the previous subsection, but now for the y and 
r variables, it is straightforward to show that the general solution to this 
equation with the initial condition AF(0, r) = has the form 

af — - r y d y' r ^ r u f 2 ^ f) , m 



2ni J 2ni uof — 4 a s 

where F 2 is the Laplace transform in y and r of the function F 2 (y,r). The 
contours for the / and u integrals lie to the right of all singularities in these 
variables. We now need a reasonable estimate for F 2 . Note that F 2 corre- 
sponds to the (Laplace tranform of) the properly normalized initial gluon 
distribution at low virtuality and small y. A rough estimate can be made 
usine the methods of ||. 

A 2 

F 2 = (48) 

[u - 2u ][f - 2 7 (u = lu )} ' V ; 

where Aq is the normalization factor for the gluon structure function, co>o is 
defined in eq. (1) and 7(0;) = — is the anomalous dimension for the leading 
twist operator. Here we consider the EKL solution to parametrize the data 
over a wide kinematic region, including GLR nonlinear corrections. Shortly 
we will discuss the case where both the GLR and multigluon corrections are 
considered small. 

All contours in (|47|) are to the right of all singularities in u and /. The 
integrand has two poles in u, one corresponding to the initial condition (2uq) 
and one from the equation (4a s //). We will demonstrate in the next subsec- 
tion that the former is dominant for the choice c^o — 0.5, and that restricting 
ourselves to its contribution is very good approximation. Under this assump- 
tion we perform the uj and / integrals, and obtain 

Ag = -4ra s 5 2 A 2 G e 2v [ V dy'e 2uov ' + 2 ^ = uo)r . (49) 

Jo 

The factor r in front arises from the double pole in the / variable. The 
answer clearly satisfies the boundary conditions. We derive further 

Aff(y,r,7y) _ 2a 8 r 2 _ F 2 (y,r) - F 2 (y = 0,r) 
g(y,r,rj) ~ ujq F(y,r) 
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This implies 



Ax B G(x B , Q ) 2a s r 2 1 

: L — () A . f ' ' ' 1 

x B G(x B ,(5 2 ) cj Q 2 R 2 n 



S 2 -t^tXbG(x b ,Q 2 ) , (51) 



if the value of the structure function is large enough in the region of low x B . 
Recall that the correct definition of rj is e v = 1/(Q 2 R%). Substituting in 
eq. ( [5lD the value of gluon structure function from HERA data |§ at Q 2 = 
15 GeV 2 and x B = 10" 4 (x B G ~ 30) and a typical value of R 2 N = 5GeV~ 2 
we obtain (with a s = 0.25 ) 

Ax B G(x B ,Q 2 ) 2 3 
— __ ^ OAb ~ 6 • 10 . 

x B G(x B ,Q 2 ) 

It is more instructive to compare the above correction to the gluon struc- 
ture function with the one due to the GLR equation. In this case we consider 
both the correction to the GLAP equation due to the GLR shadowing and 
due to multigluon correlations as small. Thus we try to find the solution to 
the GLR equation in the form: 

F{y,r) = F GLAP + AF GLR . 

For AF GLR we can write the equation: 

dydr ~ asAFGLR ~ QW N FGLAP (52) 

The solution is 

AFglr = _ * CuffM.. (*L ^ + (53) 



Q 2 R% Jo J 2-7T2 J 2m ujf — as 

where F GLAP is the Laplace image of function F GLAP (y,r). Again using (pf) 
we get 

AFqlr = - Q 2 ^fj 2 A 2 G [ V dy'e 2 ^' [e 2 ^=^ - e^=^} (54) 
6asQ z R N Jo 

Assuming that the second term in the above equation is much smaller than 
the first we have 

AFglr = ~ 3Q 2 R 2 N N C Fglap (55) 
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Finally, we can get for the ratio 



Ax B G(x B , Q 2 ) _ S& c HQ 2 IQl)p 



A(x b G(x b ,Q 2 ))glr tt 2 u;o7 

which gives a value of the order of 0.04 if ujq ~ 0.5. We will return to a 
discussion of the corrections to the GLR equation in the next section where 
we will consider the general solution to the new equation at fixed as- 



The above estimates seem in contradiction with the estimates in [fL4 
where a large contribution from multigluon correlations was found. The 
method we use here is however quite diffent from the one in [14]. There 



the effect of including the new (pole) singularity in / (the Laplace conjugate 
variable to r) that results from the resummation of bubbles associated with 
the 4-Pomeron coupling was contrasted with the contribution from the two 
Pomeron cut at the level of the Green functions. Although the location of the 
singularities is quite close, their nature and residues are very different. This 
led to a large ratio of the contributions of these singularities to the Green 
functions. In the present case we include in our estimates the two-gluon 
source, i.e. the initial gluon distribution, and renormalize it in both cases to 
the same physical initial condition (here the EKL ansatz), thus absorbing the 
residues. As was remarked in |TJ]], one can absorb the residues alternatively 
in Rn- Further, in closing the contours involved in performing the inverse 
Laplace transforms we closed on the singularities of the initial condition (be- 
ing the rightmost singularities), and not on the propagator poles. Therefore 
our renormalized R 2 N has no extra dependence on ln(l/x). We believe that 
our method is in the above sense more physical. 

5.4 Numerical Estimates 

In this subsection we will solve eq. (flEp numerically. The method we use 
goes as follows. We first perform a Laplace transform with respect to y on 
(PU|), which leads to 

dAF(u,r) 4a s Aa s 5 2 2 

a = AF(uj,r) F (uj,r) . (57) 

or uj uj 

Using various ansdtze for F(y, r) we solve this equation using Numerov's 
method, and perform finally the inverse Laplace transform with respect to 
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uo For the various ansatze we fit a parametrized form to F 2 (y, r) of which 
the Laplace transform can be taken analytically. 

We will investigate three cases. The first is the EKL ansatz from the 
previous subsection 

F(y,r) = A G e^y + f° r , (58) 

where we take uq = 0.5, fo = ds/^o, and as — 0.25. To correspond with 
the numbers given in the previous subsection we put A g ~ 0.07. We will use 
this case to check the accuracy of the estimate given earlier, and drop for 
this check the requirement AF(y = 0, r) = 0. Let us denote the numerical 
answer by AF EKL (num). Then for the pole answer we find 

AF EKL (pole) = r e 2 ^+ 2 ^ , (59) 

and for the leading term of the saddle point contribution 



AF EKL (saddle) ^. ^ 



./'; 



3 



loq \ 16Tra s y (fs - 27, 



oj 2 



where fs = \J^aa,y/r. Note that both contributions vanish as r — > 0. In Ta- 
ble 1 we list these three contributions for various x at a fixed Q 2 = 15 GeV/c 2 . 
It is clear from this table that AF EKL (pole) is a very good approximation to 
AF EKL (num), in fact better than one might expect from AF EKL (saddle). 

Next we treat a more realistic case. We now use for F(y,r) alternatively 
the MRSD0' and MRSD-' [TE\ gluon distribution functions, which are, at the 
starting scale Qq, constant as function of x, and behave as x~ ' 5 respectively. 
They are both parametrized by 

xG(x,Q 2 ) = A G x Xg (l -x) Vb (1 + i g x) , (61) 



with the coefficients AG,X g ,r] g and 7 5 given at Q = Q in [T5| . We kept 
this parametrization up until Q 2 = 15 GeV/c 2 , but refitted the coefficients 
for every step in the Numerov procedure. The Laplace transform of this 
parametrization is easily determined. Following the methods described in 
the above we determine AF(y,r). We checked that AF vanishes for small 
y. Furthermore eq. ( p7|) was solved under the condition AF(y,r = 0) = 0. 

5 We thank Keith Ellis for providing subroutines that perform this last step. 



18 



Similarly to the previous subsection we determined at Q 2 = 15 GeV/c 2 
for values of x B from 0.01 down to the LEP Cg> LHC value of 10~ 5 the ratios 



A(x B G(x Bl Q 2 )) 1 A%r) 
x B G(x B ,Q 2 ) Q 2 R 2 N x B G(x B ,Q 2 ) 1 } 

and 

A(x B G(x B ,Q 2 )) _ AF(y,r) 
A(x B G(x B ,Q 2 )) GLR ah(x B G(x B ,Q 2 )y 1 ' 

These ratios are given in Table 2. We infer from this table that corrections to 
the gluon structure function from multigluon correlations beyond the next- 
to-leading twist are small, at most 5% at small x. We see that the MRSD-' 
distribution leads to larger corrections than the MRSDO' one. As a fraction 
of the GLR correction both ansdtze are small, about 10% for the MRSD-' 
case and 10-30% for the MRSDO' case. 

Thus we confirm here the rough estimates from the previous section. 
However we note that estimates of the gluon correlation radius range 
from 1 fm to 0.3 fm. The numbers in table 2 are accordingly easily adjusted. 



6 The General Solution (for Fixed ag) 



In this section we will discuss the general solution to eq. (|17|) and its conse- 
quences for the gluon structure function. 

We start by noting that the equation (|T7|) can be written, at fixed as, as 



+ = C 2 „ . *» G C»>(z B , r + ,„) - (64) 

oyo(r + rjo) 

n - a 2 s ^e- {r+Vo) x^ +1) G {n+l \y , r + Vo ), 

where we used the fact that that x B G^ only depends on \n(Q 2 /Ql) with 
an arbitrary Q . The above form of the equation reflects the choice t] = 
— hiQlR 2 N . We focus now on the hypersurface r/ = rj. The evolution equa- 
tion (see eq. (|20D ) can then be written in the form: 



d 2 g(y,£,v) - d 2 g a s 5 2 d 4 g d 2 g 2 , dg 

dydi = as ^f + -^r ( w ~ W ] ~ a * 1 e % ~ 9) ' (65) 
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where £ = r + rj. The advantage of this form is that all explicit rj dependence 
such as exp(— rj) has been removed. 

We can find g(y,£ = r + r],rj) using a double Laplace transform with 
respect to y and r], namely 

g(y,Z,v) = j^^ eU;y + pr> 9^,P,0- (66) 

The function g((jJ,£,p) obeys the equation 

u — = {a s p + -g-(P -P ) - «s7(P - l)e 5 }g{u,£,P) (67) 

The solution to eq. ( pTf ) is: 



j (27U)" 1 

(68) 

The function g(u,p) must be determined from the initial condition, eq. fl22[) , 
for |r/| <C r, r ^> 1 and y = where the solution looks as follows: 

9<V = 0,r,v) = f ^ 2 g(uJ,p) ■ e^ e ^H P ) + 4^( P -i)(e-^-i K (6Q) 

J [Z7tl) A 



where 

uj 3a; 



0(p) = + ^l p 2 {p+1){p _ i: 



We now assert that 

g(u,p) = r(-/(p)) • ^ e^ 2 (70) 

dp 

satisfies the initial condition of eq. fl22|) with gLLA(xB, Qq) = 5(y) at r = 
ln(Q 2 /Qo) = 0- Here /(p) = p + (f>(p) and T(—f) is the Euler gamma 
function. We will prove this shortly. 

Thus, finally, the solution to eq. (|65| ) looks as follows: 



= /^r(-/o»)).^ . e w,, { + *^«- P r + ., (71) 

(27^)^ ap 
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or changing the integration from p to /, 

9(u,P,0 = j^n-f) ■ e /^*( P (/M)— P (/)^^, (72) 

where p(f) is determined by 

/ = P(f) + <P(p(f)) ■ (73) 

The contour of integration over / is defined in a such a way that all sin- 
gularities in r(— /) except the one at / = are located to the right of the 
contour. 



We can now verify the claim made in eq. (|70|). At y = y , \r]\ <C r, 
( ~ r > 1 we have 

9(Vo, r,v)^J^ 2 e^T(-f) e^ (/m/ - (74) 

We now close the / contour in the right half plane. Neglecting the 5 2 term 
we have 4>{p{n)) ~ ot$n 2 fu, and we get indeed 

g(y , r, V )^Yl / fre*™* 2 ^ (75) 

n=l U - J 11X1 

The structure function can be found from eq. fl72f) putting £ = (see eq. (|19|)). 
We obtain 

0(v,* = O,r) = /^r^.c^-Kflr^v. (76) 

Eq. (^|) implies that p(f) -> (f^)i Therefore we can close the contour 
in / over the singularities of r(— /). Next we must integrate over uj. We 
evaluate this integral using the method of steepest descent and the large / 
approximation for p(f). Neglecting terms proportional to as in the exponent 
we find the saddlepoint uj = {por/AyY^ where po = (?>n/ 6>s5 2 ) l l A . Then 



9fat=0,r) = £ L 4 L C7(n,y)e- 4 " 

n=l n - 



where C(n,y) is a pre-exponential, smooth factor. The above series clearly 
converges and the typical value of n in this series is of the order of unity. We 
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thus confirm the estimates from the previous section, and see that we could 
trust our calculations of the anomalous dimension 72 n . We note that this 
series has an infinite radius of convergence. Thus for the simplified model we 
consider and for the case of an eikonal initial condition we conclude that there 
is no analogy to a renormalon in the Wilson Operator Product Expansion. 

Now we wish to consider the solution near / — ► 1, to understand under 
which circumstances it suffices to take only this singularity into account. 
Note that this corresponds to the leading twist case, with GLAP evolution. 
In this region we can rewrite the solution in the following form. Substituting 
/ = 1 + t we obtain: 

where 

W = toy — r H r + tr\ — + —2— — 1} (78) 

oj oj ojr 

From this form of the exponent we see that for uj smaller than oj cr where 

"cr = 2a s (l + f )•{!+ ( fll 2l ^ + ( 79 ) 

3 2a s (l + /3) r uj%.r 

we cannot close the contour in t over singularities at positive t. (Here we 
introduce the parameter oj® r in order to separate the effects of 7 and 5.) For 
such oj one would need all singularities in /. To understand what happens 
in this region let us expand ^ by writing oj = oj cr + A: 

* = oj cr y - r H -r + A(y - ) (80) 

0J cr 0J cr 0J cr 

Integration over A [] gives rise to the delta function 8[-^-(t — t ) ] where 

to = —^—{oj 2 cr y - a s r] (81) 

0J cr r 

Carrying out the integration gives 

0(y f £ = O,r) = ^.±.e^y-r + ^r (g2) 



6 Alternatively one may expand to 0(A 2 ) and use steepest descent. 
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Let us define the "critical line" |4[] by 



The gluon structure function is on this critical line 

x B G(x B , g 2 ) = R% g 2 ^ £ + J" 1 (84) 

2 asr 3 47ro;s r 



For S = this becomes 



2 /V 2 

x B G(x B ,Q 2 )=I$ f Q 2 f-. (85) 



Note that this is similar to the solution of the GLR equation with running 
coupling H], but not quite the same. 

Thus, the structure of the solution to the new evolution equation with 
the initial condition eq. (|22| ) looks as follows. In the kinematic region to the 
right of the critical line we can in fact solve the linear GLAP equation, but 
with the new initial condition eq. fl8"5] ) on the critical line. To the left of the 
critical line we need the solution to the full equation. Note that if we would 
change the initial condition of the full evolution equation the value of the 
structure function on the critical line would also change. 

We have thus achieved a new understanding of the role of the initial 
condition in the problem. In particular, we conclude that the solution on the 
critical line depends only on the initial condition in the region of / — > 1, i.e. 
it depends only on the initial condition for the GLAP equation. Implicitly we 
have used here the assumption, expressed in our initial condition (p2|) , that 
the multigluon correlations are sufficiently small at large x. We recall that the 
original derivation of the GLR equation was based on this same assumption. 
We are not restricted to such an assumption for our equation. Clearly, if 
there were strong correlations between gluons at large x it would change the 
explicit form of the solution of our evolution equation. Nevetheless, the line 
of reasoning followed in this section would continue to hold. 

For 5 we have two different situations. In the first, for 

— < 



3 2a$ r 
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the solution on the critical line is the same as in eq. (B5|). I.e. the only change 
that occurs in the solution of the case 5 = is that there is a new equation 
for the critical line, eq. fl8"3"|). Note that the HERA experiments correspond 
to this situation. 

At very large values of r (Q 2 ^> Qq) , when 



3 2asr' 
the solution on the critical line looks as follows: 

x B G(x B ,Q 2 ) = R 2 N Q 2 p^- (86) 

Z7T0 Z r 

In this case we must solve the GLAP equation using eq. ( |S6| ) as the boundary 
condition. 



7 Conclusions. 

The main result of the paper is the new evolution equation ( POD . It allows us 
to penetrate deeper into the region of high density QCD because we incor- 
porate multigluon correlations into the evolution, and to answer questions 
which could not be answered before. For example one could now investigate 
the question of how well the Glauber theory for shadowing corrections in 
deep-inelastic scattering with a heavy nucleus works. 

This present equation solves two theoretical problems which arise in the 
region of high parton density: (i) It takes induced multigluon correlations into 
account, which originate from parton-parton (mainly gluon-gluon) interac- 
tions at high enegy and can be calculated in the framework of perturbative 
QCD, and (ii) it allows for an arbitrary initial gluon distribution, which is 
nonperturbative in nature. 

We have found the general solution to the new equation for the case 
of an eikonal initial condition and fixed a s . We found no evidence for a 
"renormalon" in the twist expansion. 

Our numerical estimates show that the effect of multigluon correlations is 
rather small in the accessible region of energy. We have seen evidence for this 
by using approximate methods and the general solution to the new equation. 
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We have shown that the general solution confirms the strategy developed 
for the GLR equation: we have calculated the new critical line for the gener- 
alized equation and shown that to the right of this critical line we can solve 
the linear GLAP equation with a new boundary condition on this line. We 
found this boundary condition taking into account the multigluon correla- 
tions. This approach, developed in this paper simplifies also the solution to 
the GLR equation and allows us to understand how solutions to the GLR 
equation depend on the initial conditions. This is essentially a consequence 



of the linearization of the GLR equation in ^0 



We have not discussed here the behaviour of the solution in the region 
to the left of the critical line, where multigluon correlations should come 
more forcefully into play. We plan to do this in later publication. We hope 
that the solution in the latter kinematic region will have a significant impact 
on understanding the scale of the shadowing correction and the importance 
of multigluon correlations in the so-called Regge domain. This must be 
understood in order to provide a matching between soft and hard processes. 
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X B 


AF KKL (num) 


AF EKL (poie) 


AF KKL (saddle) 


w~ 2 

lO" 4 


-3.26- 1(T 2 
-3.44 


-3.41 • 10~ 2 
-3.41 


8- 10~ 3 
2 • 1(T 2 



Table 1. Comparison of various contributions in the EKL approximation. 
Here Q 2 = 15 GeV 2 /c 2 , a s = 0.25 and R 2 N = 5 GeV" 2 . 



Ansatz 


x B 


A(x B G)/x B G 


A(x B G)/A(x B G) GLR 


MRSDO' 


lO" 2 
10~ 4 
10" 5 


-1.4 • 10" 3 
-8.9 • lO" 3 
-14.7- 10" 3 


0.11 
0.29 
0.31 


MRSD-' 


lO" 2 
10" 4 
10" 5 


-1.2 • 10~ 3 
-13.5 • lO" 3 
-40.5 • 10~ 3 


0.084 
0.1 
0.09 



Table 2. Correction to the gluon distribution function for two different 
anstitze. Here Q 2 = 15 GeV 2 /c 2 , a s ~ 0.21 and R 2 N = 5 GeV" 2 . 

FIGURE CAPTIONS 

Figure 1. 

Rescattering of Pomerons in the t-channel. 
Figure 2. 

Pictorial representation of the generalized evolution equation. 

Figure 3. 

"Fan" diagram. 

Figure 4. 

Production of three gluon shadows in a parton cascade. 



Figure 5. 



Example of type of multigluon interactions that the generalized evolution 
equation takes into account. 
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